package org.cbio.mutation_to_expression;

import org.biopax.paxtools.pattern.miner.SIFEnum;
import org.cbio.causality.analysis.*;
import org.cbio.causality.binintanalysis.Dataset1;
import org.cbio.causality.data.portal.BroadAccessor;
import org.cbio.causality.data.portal.CBioPortalAccessor;
import org.cbio.causality.data.portal.ExpDataManager;
import org.cbio.causality.model.Alteration;
import org.cbio.causality.model.AlterationPack;
import org.cbio.causality.model.Change;
import org.cbio.causality.network.*;
import org.cbio.causality.util.*;

import java.awt.*;
import java.io.*;
import java.util.*;
import java.util.List;

/**
 * This class searches for links from mutated genes in a genomic dataset to the regulated downstream
 * expressions.
 *
 * @author Ozgun Babur
 */
public class MutationToExpression
{
	/**
	 * False discovery rate threshold for results.
	 */
	private static double FDR_THR = 0.05;

	/**
	 * Threshold for log standard dev of target gene expression.
	 */
	private static double LOG_SD_THR = 4;

	/**
	 * Genomic dataset to use.
	 */
	private static final Dataset1 DATASET = Dataset1.BRCA;

	/**
	 * If a comparison has less than 2 values at one side, then we cannot use t-test, so we do not
	 * make this comparison.
	 */
	private static final int MIN_GROUP_SIZE = 2;

	/**
	 * Length limit for the paths used for relating mutations to expression changes. This is the
	 * total length limit. Since the last edge is a transcriptional regulation, the length limit for
	 * the post-translational signaling path is one less than this.
	 */
	private static int PATH_LENGTH_LIMIT = 3;

	/**
	 * Used for ignoring regulators that are mutated less than a threshold ratio.
	 */
	private static double MIN_MUT_RATIO = 0.01;

	/**
	 * If true, then hyper-mutated and hypo-mutated samples are exluded from analysis.
	 */
	private static boolean EXCLUDE_HYPERMUTATED = false;

	/**
	 * Whether to include the merged network in network comparisons.
	 */
	private static boolean CONSIDER_MERGED_NETWORK = false;

	/**
	 * Whether to limit the tested mutated genes to the MutSig significant genes.
	 */
	private static boolean LIMIT_TO_MUTSIG = false;

	/**
	 * Use this set to limit the analysis to certain target genes.
	 */
	private static final Set<String> focusExp = null;
//	private static final Set<String> focusExp = new HashSet<String>(Arrays.asList("CDKN2A"));

	/**
	 * Use this set to limit the analysis to certain regulator genes.
	 */
	private static final Set<String> focusMut = null;

	/**
	 * Directory to generate result details.
	 */
	private static final String dir = "results-details/";

	/**
	 * Post-translational signaling network.
	 */
	private Graph travSt;

	/**
	 * Transcriptional regulation network.
	 */
	private Graph travExp;

	/**
	 * Accessor to discrete cBioPortal data.
	 */
	private static CBioPortalAccessor portalAcc;

	/**
	 * Accessor to cBioPortal gene expression values.
	 */
	private static ExpDataManager expMan;

	/**
	 * Significantly mutated genes identified by MutSig.
	 */
	private static Set<String> mutsig;

	/**
	 * Used for caching subsets of samples that marks where a gene is not mutated and mutated.
	 */
	private static Map<String, boolean[][]> diffMutMap;

	/**
	 * Used for caching subsets of samples that marks where copy number is not altered.
	 */
	private static Map<String, boolean[]> cnNotChangedMap;

	/**
	 * Used for caching pairwise mut-to-targ p-values.
	 */
	private static Map<String, Map<String, Double>> pValsMap;

	/**
	 * Marks samples that are hyper-mutated or hypo-mutated.
	 */
	private static boolean[] outlier;

	public static void main(String[] args) throws IOException
	{
		configureWithProgramParameters(args);

		// post translational signaling network
		Graph graphSt;

		// transcriptional regulation network
		Graph graphExp;

		MutationToExpression mte;

		// Pathway Commons

		graphSt = PathwayCommons.getGraph(SIFEnum.CONTROLS_STATE_CHANGE_OF);
		graphExp = PathwayCommons.getGraph(SIFEnum.CONTROLS_EXPRESSION_OF);
		mte = new MutationToExpression(graphSt, graphExp);
		Map<String, Double> resPC = mte.calcPvalDistribution(true);

		// SPIKE

		graphSt = SPIKE.getGraphPostTl();
		graphExp = SPIKE.getGraphTR();
		mte = new MutationToExpression(graphSt, graphExp);
		Map<String, Double> resSPIKE = mte.calcPvalDistribution(false);

		// SignaLink

		graphSt = SignaLink.getGraphPostTl();
		graphExp = SignaLink.getGraphTR();
		mte = new MutationToExpression(graphSt, graphExp);
		Map<String, Double> resSignaLink = mte.calcPvalDistribution(false);

		Map[] pvals;
		String[] names;

		if (CONSIDER_MERGED_NETWORK)
		{
			graphSt.merge(SPIKE.getGraphPostTl());
			graphSt.merge(PathwayCommons.getGraph(SIFEnum.CONTROLS_STATE_CHANGE_OF));
			graphExp.merge(SPIKE.getGraphTR());
			graphExp.merge(PathwayCommons.getGraph(SIFEnum.CONTROLS_EXPRESSION_OF));
			mte = new MutationToExpression(graphSt, graphExp);
			Map<String, Double> resMerge = mte.calcPvalDistribution(false);
			pvals = new Map[]{resPC, resSPIKE, resSignaLink, resMerge};
			names = new String[]{"Pathway Commons", "SPIKE", "SignaLink", "Merged"};
		}
		else
		{
			pvals = new Map[]{resPC, resSPIKE, resSignaLink};
			names = new String[]{"Pathway Commons", "SPIKE", "SignaLink"};
		}

//		printDistributions(pvals);
//		System.out.println();
		printOverlaps(pvals, names);
		System.out.println();
		printFDRSizes(pvals, names);
	}

	private static void printDistributions(Map<String, Double>[] pvals)
	{
		Histogram[] h = new Histogram[pvals.length];
		for (int i = 0; i < h.length; i++)
		{
			h[i] = loadIntoHistogram(pvals[i]);
		}
		Histogram.printAll(h);
	}

	private static Histogram loadIntoHistogram(Map<String, Double> pvals)
	{
		Histogram h = new Histogram(0.05);
		h.setBorderAtZero(true);
		for (Double v : pvals.values())
		{
			h.count(v);
		}
		return h;
	}

	private static void printOverlaps(Map<String, Double>[] pvals, String[] names)
	{
		List<String>[] select = new List[pvals.length];

		for (int i = 0; i < pvals.length; i++)
		{
			select[i] = FDR.selectUsingPolyCurve(pvals[i], FDR_THR);
		}

		System.out.println();
		CollectionUtil.printNameMapping(names);
		System.out.println();
		CollectionUtil.printVennCounts(select);
	}

	private static void printFDRSizes(Map<String, Double>[] pvals, String[] names)
	{
		double max = 0.2;
		double min = 0.01;
		double[] fdrs = new double[(int) (max / min)];
		for (int i = 0; i < fdrs.length; i++)
		{
			fdrs[i] = min * (i + 1);
		}

		int[][] sizes = new int[pvals.length][];
		for (int i = 0; i < pvals.length; i++)
		{
			sizes[i] = FDR.getResultSizesUsingPolyCurve(pvals[i], fdrs);
		}

		for (String name : names)
		{
			System.out.print("\t" + name);
		}
		for (int i = 0; i < fdrs.length; i++)
		{
			System.out.print("\n" + fdrs[i]);

			for (int[] size : sizes)
			{
				System.out.print("\t" + size[i]);
			}
		}
	}

	private static void configureWithProgramParameters(String[] args)
	{
		for (int i = 0; i < args.length; i++)
		{
			if (!args[i].startsWith("-")) continue;

			if (args[i].equals("-fdr") && args.length > i+1)
			{
				FDR_THR = Double.parseDouble(args[i+1]);
			}
			else if (args[i].equals("-path-length-limit") && args.length > i+1)
			{
				PATH_LENGTH_LIMIT = Integer.parseInt(args[i + 1]);
			}
			else if (args[i].equals("-exclude-hypermutated") && args.length > i+1)
			{
				EXCLUDE_HYPERMUTATED = args[i+1].equals("on");
			}
			else if (args[i].equals("-consider-merge") && args.length > i+1)
			{
				CONSIDER_MERGED_NETWORK = args[i+1].equals("on");
			}
			else if (args[i].equals("-log2sd-threshold") && args.length > i+1)
			{
				LOG_SD_THR = Double.parseDouble(args[i + 1]);
			}
			else if (args[i].equals("-limit-to-mutsig") && args.length > i+1)
			{
				LIMIT_TO_MUTSIG = args[i+1].equals("on");
			}
			else if (args[i].equals("-min-mutation-ratio") && args.length > i+1)
			{
				MIN_MUT_RATIO = Double.parseDouble(args[i+1]);
			}
		}

		System.out.println("-fdr " + FDR_THR);
		System.out.println("-path-length-limit " + PATH_LENGTH_LIMIT);
		System.out.println("-exclude-hypermutated " + (EXCLUDE_HYPERMUTATED ? "on" : "off"));
		System.out.println("-consider-merge " + (CONSIDER_MERGED_NETWORK ? "on" : "off"));
		System.out.println("-log2sd-threshold " + LOG_SD_THR);
		System.out.println("-limit-to-mutsig " + (LIMIT_TO_MUTSIG ? "on" : "off"));
		System.out.println("-min-mutation-ratio " + MIN_MUT_RATIO);
		System.out.println();
	}

	public MutationToExpression(Graph graphSt, Graph graphExp)
		throws IOException
	{
		this.travSt = graphSt;
		this.travExp = graphExp;

		System.out.println("\nPost-translation signaling network: " + travSt.getName());
		System.out.println("Transcriptional regulation network: " + travExp.getName());

		initDataAccessors();
		initOutliers();
	}

	private void initDataAccessors() throws IOException
	{
		if (portalAcc == null) portalAcc = new CBioPortalAccessor(DATASET.mutCnCallExpZ);

		if (expMan == null)
		{
			expMan = new ExpDataManager(portalAcc.getGeneticProfileById(DATASET.exp.profileID[0]),
				portalAcc.getCaseListById(DATASET.exp.caseListID));
			expMan.setTakeLog(true);
		}
		if (LIMIT_TO_MUTSIG && mutsig == null) mutsig = BroadAccessor.getMutsigGenes(DATASET.name(), 0.05);
		if (pValsMap == null) pValsMap = new HashMap<String, Map<String, Double>>();
		if (diffMutMap == null) diffMutMap = new HashMap<String, boolean[][]>();
		if (cnNotChangedMap == null) cnNotChangedMap = new HashMap<String, boolean[]>();
	}

	private void initOutliers()
	{
		if (EXCLUDE_HYPERMUTATED && outlier == null)
		{
			Set<String> symbols = PathwayCommons.getGraph(SIFEnum.CONTROLS_STATE_CHANGE_OF).getSymbols();
			symbols.addAll(SPIKE.getGraphPostTl().getSymbols());
			symbols.addAll(SignaLink.getGraphPostTl().getSymbols());
			symbols.addAll(PathwayCommons.getGraph(SIFEnum.CONTROLS_EXPRESSION_OF).getOneSideSymbols(true));
			symbols.addAll(SPIKE.getGraphTR().getOneSideSymbols(true));
			symbols.addAll(SignaLink.getGraphTR().getOneSideSymbols(true));

			Set<AlterationPack> packs = new HashSet<AlterationPack>();

			for (String gene : symbols)
			{
				AlterationPack pack = portalAcc.getAlterations(gene);
				if (pack != null) packs.add(pack);
			}

			outlier = Outlier.getOutlierAltered(packs, Alteration.MUTATION);
		}
		else if (outlier == null)
		{
			outlier = new boolean[portalAcc.getCurrentCaseList().getCases().length];
		}
	}

	/**
	 * Gets the significant regulator target pairs as strings, like "Gene1 - Gene2".
	 */
	public Map<String, Double> calcPvalDistribution(boolean generateGraphs) throws IOException
	{
		Map<String, Set<String>> targUp2mut = new HashMap<String, Set<String>>();
		Map<String, Set<String>> targDw2mut = new HashMap<String, Set<String>>();
		Set<String> muts = new HashSet<String>();
		Set<String> tars = new HashSet<String>();

		List<String[]> pairs = getPairsToTest();
		System.out.println("test size = " + pairs.size());

		Map<String, Double> pvals = new HashMap<String, Double>();

//		Progress prg = new Progress(pairs.size());
		for (String[] pair : pairs)
		{
//			prg.tick();
			double pv = calcDiffPval(pair[0], pair[1]);
			if (!Double.isNaN(pv)) pvals.put(pair[0] + " " + pair[1], pv);
		}

//		List<String> select = selectOverThr(pvals, 0.01);
		List<String> select = FDR.selectUsingPolyCurve(pvals, FDR_THR);

		for (String key : select)
		{
			String[] name = key.split(" ");

			String mut = name[0];
			String tar = name[1];

			if (!targUp2mut.containsKey(tar)) targUp2mut.put(tar, new HashSet<String>());
			if (!targDw2mut.containsKey(tar)) targDw2mut.put(tar, new HashSet<String>());

			muts.add(mut);
			tars.add(tar);

			boolean up = calcMeanChange(tar, Collections.singleton(mut)) > 0;
			if (up) targUp2mut.get(tar).add(mut);
			else targDw2mut.get(tar).add(mut);
		}

		System.out.println("result size = " + select.size());
		System.out.println("muts = " + muts.size());
		System.out.println("targs = " + tars.size());

		System.out.println("CDKN2A upstr: " + targDw2mut.get("CDKN2A"));

		if (generateGraphs) generateResultDetails(targUp2mut, targDw2mut);

		return pvals;
	}

	private List<String> selectOverThr(Map<String, Double> pvals, double thr)
	{
		List<String> select = new ArrayList<String>();
		for (String names : pvals.keySet())
		{
			double pv = pvals.get(names);
			if (pv <= thr)
			{
				select.add(names);
			}
		}
		return select;
	}

	private List<String[]> getPairsToTest()
	{
		List<String[]> pairs = new ArrayList<String[]>();
		Set<String> approvedMuts = new HashSet<String>();

		Set<String> tars = focusExp != null ? focusExp : travExp.getOneSideSymbols(false);

		for (String tar : tars)
		{
			double[] v = expMan.get(tar);
			if (v == null) continue;
			double[] vals = new double[v.length];
			for (int i = 0; i < v.length; i++)
			{
				vals[i] = Math.pow(2, v[i]);
			}
			double sd = Math.log(Summary.stdev(vals) / Math.log(2));
			if (sd < LOG_SD_THR) continue;

			Set<String> muts = getUpstream(tar);
			if (focusMut != null) muts.retainAll(focusMut);
			else if (LIMIT_TO_MUTSIG) muts.retainAll(mutsig);

			for (String mut : muts)
			{
				if (!approvedMuts.contains(mut))
				{
					AlterationPack pack = portalAcc.getAlterations(mut);
					if (pack == null) continue;
					if (pack.get(Alteration.MUTATION) == null) continue;
					if (pack.getAlteredRatio(Alteration.MUTATION) < MIN_MUT_RATIO) continue;

					approvedMuts.add(mut);
				}
				pairs.add(new String[]{mut, tar});
			}
		}
		return pairs;
	}

	/**
	 * Gets upstream of the given target gene. Makes sure that the relating paths contain
	 * transcriptional regulation relations only at the last step.
	 */
	private Set<String> getUpstream(String targ)
	{
		Set<String> expUp = new HashSet<String>(travExp.getUpstream(targ));

		if (PATH_LENGTH_LIMIT > 1) expUp.addAll(travSt.getUpstream(expUp, PATH_LENGTH_LIMIT - 1));

		return expUp;
	}

	/**
	 * Gets the not mutated and mutated subsets of the given gene. Also makes sure that samples are
	 * not altered in copy number or expression.
	 */
	private boolean[][] getDiffMutatedSubsets(String gene)
	{
		if (diffMutMap.containsKey(gene)) return diffMutMap.get(gene);

		boolean[][] diff = null;

		AlterationPack alts = portalAcc.getAlterations(gene);
		if (!(alts == null || alts.get(Alteration.MUTATION) == null ||
			alts.getAlteredRatio(Alteration.MUTATION) < MIN_MUT_RATIO))
		{
			diff = new boolean[2][];
			diff[0] = selectSubset(Collections.singleton(gene), false);
			diff[1] = selectSubset(Collections.singleton(gene), true);
		}

		diffMutMap.put(gene, diff);
		return diff;
	}

	/**
	 * Gets the samples that the gene is either mutated or intact. In both cases it is expressed,
	 * and not copy number lost.
	 */
	private boolean[] selectSubset(Set<String> symbols, boolean mutated)
	{
		AlterationPack[] alts = new AlterationPack[symbols.size()];

		int i = 0;
		for (String symbol : symbols)
		{
			alts[i++] = portalAcc.getAlterations(symbol);
		}

		boolean[] x = new boolean[alts[0].getSize()];

		for (i = 0; i < x.length; i++)
		{
			x[i] = (mutated ? atLeastOneIsMutated(alts, i) : nothingMutated(alts, i)) &&
				copyNumberNotAltered(alts, i) && expressionNotAltered(alts, i) &&
				!outlier[i];
		}

		return x;
	}

	private boolean atLeastOneIsMutated(AlterationPack[] alts, int index)
	{
		for (AlterationPack alt : alts)
		{
			if (alt.getChange(Alteration.MUTATION, index).isAltered()) return true;
		}
		return false;
	}

	private boolean nothingMutated(AlterationPack[] alts, int index)
	{
		for (AlterationPack alt : alts)
		{
			if (alt.getChange(Alteration.MUTATION, index) != Change.NO_CHANGE) return false;
		}
		return true;
	}

	private boolean copyNumberNotAltered(AlterationPack[] alts, int index)
	{
		for (AlterationPack alt : alts)
		{
			Change ch = alt.getChange(Alteration.COPY_NUMBER, index);
			if (ch != Change.NO_CHANGE) return false;
		}
		return true;
	}

	private boolean expressionNotAltered(AlterationPack[] alts, int index)
	{
		for (AlterationPack alt : alts)
		{
			Change ch = alt.getChange(Alteration.EXPRESSION, index);
			if (ch != Change.NO_CHANGE) return false;
		}
		return true;
	}

	private double[][] getValueSubsets(String symbol, boolean[] set1, boolean[] set2)
	{
		double[] exp = expMan.get(symbol);

		if (exp == null) return null;

		boolean[] cnUnchanged = getCopyNumberUnchanged(symbol);
		if (cnUnchanged == null) return null;

		double[] vals1 = getSubset(exp, set1, cnUnchanged);
		double[] vals2 = getSubset(exp, set2, cnUnchanged);
		return new double[][]{vals1, vals2};
	}

	private double calcDiffPval(String mut, String targ)
	{
		if (pValsMap.containsKey(mut) && pValsMap.get(mut).containsKey(targ))
			return pValsMap.get(mut).get(targ);

		double pval;
		boolean[][] diff = getDiffMutatedSubsets(mut);

		if (diff == null)
		{
			pval = Double.NaN;
		}
		else
		{
			double[][] vals = getValueSubsets(targ, diff[0], diff[1]);
			if ( vals == null || vals[0].length < MIN_GROUP_SIZE || vals[1].length < MIN_GROUP_SIZE)
			{
				pval = Double.NaN;
			}
			else
			{
				pval = calcPval(vals[0], vals[1]);

			}
		}

		if (!pValsMap.containsKey(mut)) pValsMap.put(mut, new HashMap<String, Double>());
		if (!pValsMap.get(mut).containsKey(targ)) pValsMap.get(mut).put(targ, pval);

		return pval;
	}

	private double calcDiffPval(String symbol, boolean[] set1, boolean[] set2)
	{
		double[][] vals = getValueSubsets(symbol, set1, set2);
		if ( vals == null) return Double.NaN;

		if (vals[0].length < MIN_GROUP_SIZE || vals[1].length < MIN_GROUP_SIZE) return Double.NaN;

		return calcPval(vals[0], vals[1]);
	}

	private double calcMeanChange(String symbol, boolean[] set1, boolean[] set2)
	{
		double[][] vals = getValueSubsets(symbol, set1, set2);
		if ( vals == null) return Double.NaN;

		if (vals[0].length < MIN_GROUP_SIZE || vals[1].length < MIN_GROUP_SIZE) return Double.NaN;

		return Summary.calcChangeOfMean(vals[0], vals[1]);
	}

	private double calcPval(double[] v1, double[] v2)
	{
		return StudentsT.getPValOfMeanDifference(v1, v2);
	}

	private double[] getSubset(double[] vals, boolean[] inds, boolean[] cnUnchanged)
	{
		List<Double> list = new ArrayList<Double>(vals.length);

		for (int i = 0; i < vals.length; i++)
		{
			if (inds[i] && cnUnchanged[i] && !outlier[i] &&
				!Double.isNaN(vals[i])) list.add(vals[i]);
		}
		double[] sub = new double[list.size()];
		for (int i = 0; i < sub.length; i++)
		{
			sub[i] = list.get(i);
		}
		return sub;
	}

	private boolean[] getCopyNumberUnchanged(String gene)
	{
		if (cnNotChangedMap.containsKey(gene)) return cnNotChangedMap.get(gene);

		AlterationPack alts = portalAcc.getAlterations(gene);

		if (alts == null) return null;

		Change[] cnc = alts.get(Alteration.COPY_NUMBER);

		if (cnc == null) return null;

		boolean[] b = new boolean[cnc.length];
		for (int i = 0; i < b.length; i++)
		{
			b[i] = !cnc[i].isAltered() && !cnc[i].isAbsent();
		}

		cnNotChangedMap.put(gene, b);

		return b;
	}

	/**
	 * Generates influence graphs for targets with at least 4 regulators.
	 */
	private void generateResultDetails(Map<String, Set<String>> targUp2mut,
		Map<String, Set<String>> targDw2mut)
	{
		for (String targ : targUp2mut.keySet())
		{
			Set<String> muts = targUp2mut.get(targ);
			if (muts.size() > 2) generateInfluenceGraphs(targ, muts, true);
		}
		for (String targ : targDw2mut.keySet())
		{
//			if (!targ.equals("CDKN2A")) continue;
			Set<String> muts = targDw2mut.get(targ);
			if (muts.size() > 2) generateInfluenceGraphs(targ, muts, false);
		}

		System.out.println("influence graphs generated");
	}

	//--- Section: Graph generation----------------------------------------------------------------|

	/**
	 * Gets the p-val for the link between expression of the given target and mutations in any of
	 * the given regulators.
	 */
	private double calcPVal(String target, Set<String> ups)
	{
		boolean[] normal = selectSubset(ups, false);
		boolean[] mutated = selectSubset(ups, true);

		return calcDiffPval(target, normal, mutated);
	}

	/**
	 * Gets the mean difference of the target expression when the given regulators are mutated
	 * compared to not mutated.
	 */
	private double calcMeanChange(String target, Set<String> ups)
	{
		boolean[] normal = selectSubset(ups, false);
		boolean[] mutated = selectSubset(ups, true);

		return calcMeanChange(target, normal, mutated);
	}

	/**
	 * Gets the two sets of values of target expression to compare.
	 */
	private double[][] getValueSubsets(String target, Set<String> ups)
	{
		boolean[] normal = selectSubset(ups, false);
		boolean[] mutated = selectSubset(ups, true);

		return getValueSubsets(target, normal, mutated);
	}

	/**
	 * Generates an influence graph for the given target and given upstream regulators. The graphs
	 * are saved on the disk.
	 */
	public void generateInfluenceGraphs(String target, final Set<String> muts,
		final boolean upregulation)
	{
		UpstreamTree tree = new UpstreamTree(travSt, travExp, new BranchDataProvider()
		{
			final Color TARG_UP_COLOR = new Color(220, 255, 220);
			final Color TARG_DW_COLOR = new Color(255, 255, 200);

			@Override
			public Color getColor(String gene, String root)
			{
				if (gene.equals(root))
				{
					return upregulation ? TARG_UP_COLOR : TARG_DW_COLOR;
				}
				else if (muts.contains(gene))
				{
					double pval = calcPVal(root, Collections.singleton(gene));
					String color = val2Color(pval);
					String[] c = color.split(" ");
					return new Color(Integer.parseInt(c[0]), Integer.parseInt(c[1]),
						Integer.parseInt(c[2]));
				}
				else return Color.WHITE;
			}

			@Override
			public double getThickness(GeneBranch branch, String root)
			{
				Set<String> genes = branch.getAllGenes();
				genes.retainAll(muts);

				if (genes.isEmpty())
				{
					System.out.println();
				}
				assert !genes.isEmpty();
				double pval = calcPVal(root, genes);
				return Math.sqrt(-Math.log(pval));
			}
		});

		String dir = this.dir + DATASET.name() + "/";
		File d = new File(dir);
		if (!d.exists()) d.mkdirs();
		assert d.isDirectory();

		GeneBranch g = tree.getTree(target, muts, PATH_LENGTH_LIMIT);
		g.trimToMajorPaths(muts);

		RadialInfluenceTree.write(g, true, dir + target + (upregulation ? "-up" : "-dw") + ".svg");

		writeTree(g, muts, upregulation, dir);
		writeValues(target, muts, upregulation, dir);
	}

	private void writeValues(String targ, Set<String> muts, boolean upregulation, String dir)
	{
		try
		{
			BufferedWriter writer = new BufferedWriter(new FileWriter(dir + targ +
				(upregulation ? "-up" : "-dw") + "-vals.txt"));
			writer.write("normal\tmutated");
			double[][] vals = getValueSubsets(targ, muts);

			for (int i = 0; i < Math.max(vals[0].length, vals[1].length); i++)
			{
				writer.write("\n");
				if (i < vals[0].length) writer.write("" + vals[0][i]);
				writer.write("\t");
				if (i < vals[1].length) writer.write("" + vals[1][i]);
			}

			writer.close();
		}
		catch (IOException e)
		{
			e.printStackTrace();
		}
	}

	private void writeTree(GeneBranch gwu, Set<String> muts, boolean upregulation, String dir)
	{
		try
		{
			BufferedWriter writer1 = new BufferedWriter(new FileWriter(dir + gwu.gene + (upregulation ? "-up" : "-dw") + ".sif"));
			BufferedWriter writer2 = new BufferedWriter(new FileWriter(dir + gwu.gene + (upregulation ? "-up" : "-dw") + ".format"));

			writer2.write("graph\tgrouping\ton\n");
			writer2.write("node\t" + gwu.gene + "\tcolor\t" +
				(upregulation ? "220 255 220" : "255 255 200") + "\n");

			for (GeneBranch up : gwu.branches)
			{
				String edgeTag = SIFEnum.CONTROLS_EXPRESSION_OF.getTag();
				writer1.write(up.gene + "\t" + edgeTag + "\t" + gwu.gene + "\n");

				writeWeights(gwu.gene, gwu.gene, muts, up, edgeTag, writer2);

				writePart(gwu.gene, muts, up, writer1, writer2);
			}
			writer1.close();
			writer2.close();
		}
		catch (IOException e)
		{
			e.printStackTrace();
			System.exit(1);
		}
	}

	private void writePart(String origTarget, Set<String> muts, GeneBranch gwu, Writer writer1,
		Writer writer2) throws IOException
	{
		for (GeneBranch up : gwu.branches)
		{
			String edgeTag = SIFEnum.CONTROLS_STATE_CHANGE_OF.getTag();
			writer1.write(up.gene + "\t" + edgeTag + "\t" + gwu.gene + "\n");

			writeWeights(origTarget, gwu.gene, muts, up, edgeTag, writer2);

			writePart(origTarget, muts, up, writer1, writer2);
		}
	}

	private void writeWeights(String orig, String to, Set<String> muts, GeneBranch gwu,
		String edgeType, Writer writer) throws IOException
	{
		Set<String> upstr = gwu.getAllGenes();
		upstr.retainAll(muts);
		assert !upstr.isEmpty();
		double cumPval = calcPVal(orig, upstr);

		String key = gwu.gene + " " + edgeType + " " + to;
		writer.write("edge\t" + key + "\tcolor\t" + val2Color(cumPval) + "\n");
		writer.write("edge\t" + key + "\twidth\t3\n");

		if (muts.contains(gwu.gene))
		{
			double pval = calcPVal(orig, Collections.singleton(gwu.gene));
			writer.write("node\t" + gwu.gene + "\tbordercolor\t" + val2Color(pval) + "\n");
		}
		else
		{
			writer.write("node\t" + gwu.gene + "\tbordercolor\t255 255 255\n");
		}
		writer.write("node\t" + gwu.gene + "\tborderwidth\t3\n");
		writer.write("node\t" + gwu.gene + "\tcolor\t255 255 255\n");
	}

	/**
	 * Maps p-values to gray-scale colors.
	 */
	private String val2Color(double pval)
	{
		double score = Math.min(7, -Math.log10(pval));

		int v = 250 - (int) Math.round((250D / 7D) * score);

		return v + " " + v + " " + v;
	}
}
